#clear workspace
rm(list=ls())

#load packages
packs <- c("tidyverse", "ggplot2", "stargazer", "sp", "maptools", "rgeos", "rgdal","stargazer") 
lapply(packs, library, character.only = TRUE)

#load data
village_master <- read_csv("data/village_sample_initial.csv")

df <- readRDS("data/df.rds") %>% 
  as_tibble()

# Tables ------------------------------------

################################################################
################################################################
## Table 1: Model Results
################################################################
################################################################

glm_preempt_1 <- glm(depop_preemptive ~ social_cohesion_add,data=df,family=binomial())
glm_preempt_2 <- glm(depop_preemptive ~ social_cohesion_add+vil_pop_log+min_viol_depop_dist+avg_elevation+land_area_log+total_cultivatable_land+dist_intl_border,data=df,family=binomial())
glm_preempt_3 <- glm(depop_preemptive ~ social_cohesion_add*vil_pop_log+min_viol_depop_dist+avg_elevation+land_area_log+total_cultivatable_land+dist_intl_border,data=df,family=binomial())


stargazer(glm_preempt_1,
          glm_preempt_2,
          glm_preempt_3,
          header= FALSE,
          font.size = "small",
          dep.var.labels = "Preemptive Evacuation",
          omit=c("Constant","aic"),
          covariate.labels = c("Social Cohesion",
                               "Village Pop",
                               "Dist. to Viol. Depop.",
                               "Mean Elevation",
                               "Land Area",
                               "Total Cultivatable Land",
                               "Dist. to Int'l. Border",
                               "Social Cohesion * Village Pop"),
          omit.stat=c("aic","wald","lr","ll"),
          title = "Social Cohesion and Preemptive Evacuation")



# Figures -----------------------------------

################################################################
################################################################
## Figure 1: Map of Village Evacuation Outcomes
################################################################
################################################################

villages <- readOGR("data/shapefiles/villages_id_crs3857.shp")
un_border <- readOGR("data/shapefiles/un_border_crs3857.shp")
roads <- readOGR("data/shapefiles/roads_crs3857.shp")
outside_border <- readOGR("data/shapefiles/outside_border.shp")

# prepare the 3 components: coordinates, data, and proj4string
coords <- village_master[!is.na(village_master$latitude) & !is.na(village_master$longitude) & village_master$refugee_camp == 0 & village_master$bedouin == 0& village_master$pop0 == 0 & village_master$jewish_only == 0, c("longitude", "latitude")]

# filter out villages with no population (pop0 = 1), refugee = 1, bedouin = 1
data   <- dplyr::select(filter(village_master, !is.na(latitude), !is.na(longitude), refugee_camp == 0, bedouin == 0, pop0 == 0, jewish_only == 0),-latitude, -longitude)
crs    <- CRS("+init=epsg:4326")


# make the spatial points data frame object
village_master_shape <- SpatialPointsDataFrame(coords = coords,
                                               data = data, 
                                               proj4string = crs)
village_master_shape <- spTransform(village_master_shape,CRS("+init=epsg:3857"))


village_master_shape_sf <- sf::st_as_sf(village_master_shape)


ggplot() + 
  geom_polygon(data = outside_border, aes(x = long, y = lat),colour = "black",fill=NA)+
  theme_void()+
  geom_path(data = un_border, aes(x = long, y = lat,group=id),colour = "gray20",fill=NA,alpha=.5)+
  geom_path(data = roads, aes(x = long, y = lat,group=id),colour = "gray40",fill=NA,alpha=.25, lty=2)+
  geom_sf(data=village_master_shape_sf, aes(shape = factor(depop_outcome), color = factor(included)),size=1) +
  scale_shape_discrete(name="evacuation outcome")+
  scale_colour_manual(name="areas of operation", values = c("black","grey50")) +
  theme(legend.position=c(.2, .8),
        legend.title = element_text(colour="black", size=12, face="bold"),
        legend.text = element_text(colour="black", size = 10, face = "bold"))#+

ggsave(paste("figures/Figure1_Evacuation_Outcome_Map.pdf", sep = ""), width=4.83, height = 10, units = "in")


################################################################
################################################################
## Figure 3: Distribution of Village Evacuation Outcomes
################################################################
################################################################

df_fig3 <- data.frame(main_var = c("preemptive evacuation", "remain", "remain")
                      , factor_var = c("preemptive evacuation", "reactive evacuation", "no evacuation")
                      , count=c(nrow(dplyr::filter(village_master, v_survey_collected == 1, depop_outcome == "preemptive evacuation", included == "contested"))
                                , nrow(dplyr::filter(village_master, v_survey_collected == 1, depop_outcome == "reactive evacuation", included == "contested"))
                                , nrow(dplyr::filter(village_master, v_survey_collected == 1, depop_outcome == "no evacuation", included == "contested"))
                      )
                      , order = c(3, 2, 1)
                      , order2 = c(3, 2, 1)
)


ggplot(data = df_fig3, 
       aes(x = reorder(main_var, order), y = count, fill = reorder(factor_var, order2))) + 
  geom_bar(stat="identity", position="stack") +
  labs(x = "evacuation outcome" , y = "number of villages") +
  scale_fill_grey() + 
  guides(fill=guide_legend(title="evacuation outcome")) +
  theme_light() +
  theme(legend.position=c(.775, .825), panel.grid.major = element_blank(), panel.grid.minor = element_blank())





ggsave(paste("figures/Figure3_Evacuation_Freq.pdf", sep = ""), width=4.83, height=4.75, units = "in", pointsize=11)


################################################################
################################################################
## Figure 4: Distribution of the Social Cohesion Index
################################################################
################################################################

ggplot(data=df,aes(social_cohesion_add))+
  geom_bar()+
  xlab("social cohesion index") +
  theme_classic()+
  ylab("number of villages") + 
  #ggtitle("Distribution of Social Cohesion Index Values") +
  theme_light() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

ggsave(paste("figures/Figure4_SocialCohesion_Freq.pdf", sep = ""), width=4.83, height=4.75, units = "in", pointsize=11)


################################################################
################################################################
## Figure 5: Map of Evacuation Outcome x Social Cohesion
################################################################
################################################################
village_master_shape_plus_soc <- readRDS("data/village_master_shape_plus_soc.rds")


ggplot() + 
  geom_polygon(data = outside_border, aes(x = long, y = lat),colour = "black",fill=NA)+
  theme_void()+
  geom_path(data = un_border, aes(x = long, y = lat,group=id),colour = "gray20",fill=NA,alpha=.5)+
  geom_path(data = roads, aes(x = long, y = lat,group=id),colour = "gray40",fill=NA,alpha=.25, lty=2)+
  geom_point(data=village_master_shape_plus_soc, aes(x = longitude, y = latitude, color = social_cohesion_add, shape=factor(depop_outcome,labels=c("preemptive evacuation","no evacuation","reactive evacuation"))))+
  scale_shape_discrete(name="evacuation outcome")+
  scale_colour_gradient(name="social cohesion index", low = "gray", high = "black")+
  #scale_colour_discrete(name="social cohesion index", c("gray25", "gray50", "gray75", "black"))+
  theme(legend.position = c(0.2, .8)) #+ 

ggsave(paste("figures/Figure5_Evacuation_and_SocialCohesion_Map.pdf", sep = ""), width=4.83, height=10, units = "in", pointsize=11)

################################################################
################################################################
## Figure 6: Distribution of Village Evacuation Outcomes x Social Cohesion
################################################################
################################################################


df_fig6 <- data.frame(main_var = c("0", 
                                   "0",
                                   "1",
                                   "1",
                                   "2",
                                   "2",
                                   "3",
                                   "3")
                      ,  factor_var= c("remain",
                                       "preemptive",
                                       "remain",
                                       "preemptive",
                                       "remain",
                                       "preemptive",
                                       "remain",
                                       "preemptive")
                      , count=c(nrow(dplyr::filter(df, social_cohesion_add == 0, depop_preemptive == 0,district!="Tulkarm"))
                                , nrow(dplyr::filter(df, social_cohesion_add == 0, depop_preemptive == 1,district!="Tulkarm"))
                                , nrow(dplyr::filter(df, social_cohesion_add == 1, depop_preemptive == 0,district!="Tulkarm"))
                                , nrow(dplyr::filter(df, social_cohesion_add == 1, depop_preemptive == 1,district!="Tulkarm"))
                                , nrow(dplyr::filter(df, social_cohesion_add == 2, depop_preemptive == 0,district!="Tulkarm"))
                                , nrow(dplyr::filter(df, social_cohesion_add == 2, depop_preemptive == 1,district!="Tulkarm"))
                                , nrow(dplyr::filter(df, social_cohesion_add == 3, depop_preemptive == 0,district!="Tulkarm"))
                                , nrow(dplyr::filter(df, social_cohesion_add == 3, depop_preemptive == 1,district!="Tulkarm"))
                                
                      )
)


ggplot(data = df_fig6, 
       aes(x = main_var, y = count, fill = factor_var)) + 
  geom_bar(stat="identity", position=position_dodge()) +
  labs(x = "social cohesion index" , y = "number of villages") +
  scale_fill_grey(labels = c("preemptive evacuation", "remain")) + 
  guides(fill=guide_legend(title="evacuation outcome")) +
  theme_light() +
  theme(legend.position=c(.775, .825), panel.grid.major = element_blank(), panel.grid.minor = element_blank())


ggsave(paste("figures/Figure6_Evacuation_by_Cohesion.pdf", sep = ""), width=4.83, height=4.75, units = "in", pointsize=11)

